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ABSTRACT 

Present grids of stellar atmosphere models are the workhorses in interpreting stellar observations, 
and determining their fundamental parameters. These models rely on greatly simplified models of 
convection, however, lending less predictive power to such models of late type stars. 
We present a grid of improved and more reliable stellar atmosphere models of late type stars, based on 
deep, 3D, convective, stellar atmosphere simulations. This grid is to be used in general for interpreting 
observations, and improve stellar and asteroseismic modeling. 

We solve the Navier Stokes equations in 3D and concurrent with the radiative transfer equation, for 
a range of atmospheric parameters, covering most of stellar evolution with convection at the surface. 
We emphasize use of the best available atomic physics for quantitative predictions and comparisons 
with observations. 

We present granulation size, convective expansion of the acoustic cavity, asymptotic adiabat, as func- 
tion of atmospheric parameters. These and other results are also available in electronic form. 
Keywords: Stars: atmospheres - stars: late-type - stars: interiors - stars: fundamental parameters - 
physical data and processes: convection 



1. INTRODUCTION 

We present a homogeneous grid of three dimensional 
(3D) convective stellar atmosphere models, with applica- 
tions to a wide range of astronomical inquiries. Previews 
of the grid have previously been given by Trampedach 
(2007); Trampedach & Stein (2011) with recent astero- 
seismic applications by Mathur et al. (2011, 2012) and 
Bonaca et al. (2012). The simulations have been com- 
puted with the Stein & Nordlund (1998)-code, but with 
updated atomic physics as described in more detail in 
Sect. 2. The solar abundance determination by Asplund 
et al. (2005) was based on a solar simulation identical to 
ours in all respects, except for a slightly higher horizontal 
resolution. 

The first systematic analysis of convection in stellar 
atmospheres, through 3D simulations, was performed 
in the series of papers by Nordlund & Dravins (1990); 
Dravins & Nordlund (1990). They were interested in the 
3D structure of granulation and the effects on spectral 
lines. These simulations were computed with an anelas- 
tic precursor to the present, fully compressible Stein & 
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Nordlund (1998)-code. 

Currently an independent grid, the CIFIST grid (Lud- 
wig et al. 2009), is being assembled using the C0 5 BOLD 
code (Wedemeyer et al. 2004; Freytag et al. 2012) to 
compute 3D simulations of late-type stars, including sub- 
solar metallicities and M-type stars. 

A new grid of 3D convection simulations using the 
Stagger-code is nearing completion (Collet et al. 2011; 
Magic et al. 2013), and will encompass both sub- and 
super-solar metallicities. We are truly entering the age 
of 3D convective atmosphere simulations, and in com- 
bination with recent and soon to be realized observa- 
tional advances, we are sure to learn much about the 
inner workings of stars in the coming decades. 

This paper is an attempt at publishing results of 3D 
convection simulations in a format similar to what as- 
tronomers have been using for the last half a century: 
a grid in the atmospheric parameters effective tempera- 
ture, T e fi, surface gravity, g, and metallicity, [Fe/H]. This 
grid, however, is only for solar metallicity, [Fc/H] = 0. 

A major motivation for this effort, has been the advent 
of asteroseismology and the significant leap in stellar ob- 
servations provided by the MOST (Walker et al. 2003), 
COROT (Baglin et al. 2002) and Kepler (Borucki et al. 
2010) missions, and even by the star tracker on the failed 
infra-red satellite, WIRE (e.g., Bruntt et al. 2005). These 
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observations, with their highly stable, precision photom- 
etry and long time-lines, can reveal a wealth of detail 
about field stars that could previously only be consid- 
ered for the Sun and a handful of other stars. The sheer 
volume of Kepler observations also means that a large 
range of stellar conditions are well sampled and mean- 
ingful statistics can be developed. The much anticipated 
launch in 2013 of the Gaia mission (Perryman et al. 2001; 
Prusti 2012) will greatly improve our knowledge of stel- 
lar distances, including for most Kepler targets, provid- 
ing strong and independent constraints on absolute lumi- 
nosities and radii. These observational advances demand 
reciprocal advances on the modeling side, and the present 
grid of convection simulations is our contribution to such 
an advance. 

The present work is an alternative to previously pub- 
lished grids of ID stellar atmosphere models. Some of 
the ID atmosphere grids that have found widespread 
use have been constructed by Gustafsson et al. (2008, 
MARCS), Castelli & Kurucz (2003, ATLAS9) and 
Hauschildt et al. (1999a,b, NextGen-grid of PHOENIX 
models). These models vary in atomic and molecular 
line databases, continuum opacities, equation of state, 
whether they employ plane-parallel or spherical geom- 
etry and to what extent non-LTE effects are included. 
In all these cases radiative transfer is the focus and it 
is solved in impressive detail and with high accuracy. 
For late type stars, they also all share the same short- 
coming: The analytical formulation of convection in the 
form of the widely used mixing length theory pioneered 
by Bohm-Vitense (1958, MLT) or the more recent vari- 
ant by Canuto & Mazzitelli (1991, 1992). In contrast, the 
present atmosphere grid is based on explicitly evolving 
the Navier Stoke's equations, as detailed in Sect. 2, with 
the hydrodynamics being directly coupled with realistic 
radiative transfer (see Sect. 2.3). 

The lay-out of our grid of 37 simulations, is described 
in Sect. 2.5, where the issues of interpolating in the grid, 
and methods for averaging of the simulations are also 
discussed. General properties of convection in the simu- 
lations are discussed in Sect. 3, and contrasted with the 
analytical MLT approach. Here we also touch on the 
effect of 3D convection on the frequencies of p-modes. 

The surface manifestation of convection on the Sun, 
has been observed for more than two centuries (Her- 
schel 1801), with the term granules coined by Dawes 
(1864). Granulation on other stars has been rather elu- 
sive, though, with a tentative observation on Bctclgcuse 
(aOri) by Lim et al. (1998). Our simulations reproduce 
the quiet Sun granulation, both with respect to size, con- 
trast and general shapes, giving us confidence in our pre- 
dictions for other stellar parameters, as presented in Sect. 
4. Recently, the unprecedented time-line, stability and 
precision of the Kepler observations, has allowed direct 
observations of the effect of granulation on integrated 
star light from red giants (Mathur et al. 2011). These ob- 
servations were analyzed and compared with the present 
grid of simulations, and the agreement is generally good, 
although unresolved issues remain. A more detailed anal- 
ysis, also including main sequence (MS) stars, will be 
carried out in the near future. 

We round off with conclusions in Sect. 5. 

2. THE 3D CONVECTION SIMULATIONS 



The simulations evolve the fully compressible Navier 
Stoke's equations for mass, momentum and energy con- 
servation 

d In g _ _ 
— - — = — u ■ V m g — V • u 
at 
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where g is the density, P g the gas pressure, e the inter- 
nal energy per unit mass, u the velocity field and g is the 
gravity, assumed to be constant with depth. This version 
of the equations, differ from the normal version, by hav- 
ing been divided by g to pre-condition the equations for 
the very large density gradients in stellar surface layers. 

The Stein & Nordlund (1998)-code used here is a non- 
staggered, finite difference code. Derivatives and inter- 
polations are evaluated from cubic splines and time is ad- 
vanced with a third order leapfrog (predictor/corrector) 
scheme (Hyman 1979; Nordlund & Stein 1990). 

As all higher-order (higher than linear) order schemes 
are unstable artificial diffusion is needed for stabilizing 
the solution. We employ a quenched hyper-diffusion 
scheme which efficiently concentrates the diffusion at 
sharp changes in variables, leaving the smoother parts of 
the solution unaffected. Hyper-diffusion means a fourth- 
order term is added to the normal second-order (Lapla- 
cian) term, and the quenching consists of ensuring the 
fourth-order term everywhere has the same sign as the 
second order term, to not excite more ringing. This is ex- 
plained in more detail by Stein & Nordlund (1998). This 
does not constitute a sub-grid scale model, as we have 
not specified any preconceptions of what should happen 
below the grid-scale of our simulations. That in turn 
implies that those scales do not contribute any signifi- 
cant net fluxes, i.e., motions on that scale are considered 
isotropic turbulence. 

This diffusion by the viscous stress tensor a, gives rise 
to the dissipation, <5 v i sc = ^2nj &ijdui/ drj, in Eq. (1). 
The radiative heating, Q ra( j, is dealt with in more detail 
in Sect. 2.3. 

Each of the simulations has been performed on a 
150 x 150 x 82 rectangular grid with equidistant and cyclic 
horizontal grid and a vertical (smooth) grid which is opti- 
mized to capture the large temperature gradients in the 
photosphere. The simulations cover about 13 pressure 
scale-heights vertically, with half of those below and half 
above the photosphere. On a Rosseland optical depth 
scale, this corresponds to log 10 tr OS s = [—4.5; 5-8]. This 
is the region we will refer to as stellar surface layers. The 
effects of numerical resolution on stratification and veloc- 
ity fields have been explored by Asplund et al. (2000a), 
on which we based our choice of resolution. 

Both the top and bottom boundaries are penetrable 
with outflows leaving undisturbed. To approximate the 
effect of the convection zone below the simulation do- 
main, inflows at the bottom boundary are evolved to- 
wards a constant entropy, hydrostatic equilibrium and a 
vertical velocity that balances the mass flux of the un- 
altered outflows. The inflows are evolved towards this 
ideal solution, on time-scales of half the minimum flow- 
crossing time at the bottom, which means the inflow pro- 



Grid of 3D Atmosphere Models 



3 



files are not entirely flat. The entropy is in practice set 
by values of density and internal energy. The total flux 
of the simulation, crT 4 ff , is solely due to the supplied 
entropy of the inflows — no artificial fluxes are imposed 
anywhere. Effects of the bottom boundary are imper- 
ceptible in temperature, pressure and density, but recog- 
nized in entropy, superadiabatic gradient and velocities, 
for which they affect the bottom 5-10 grid-points, cover- 
ing 0.6-1.0 pressure scale-heights. The top boundary is a 
fiducial layer, located three times further away from the 
first grid point in, than the next. This together with the 
exponentially declining density, ensures only little mass 
is present there, and minimizes its influence on the sim- 
ulation. Velocities and energies from the first grid point 
are simply copied to the fiducial layer, and the density is 
extrapolated hydrostatically. 

2.1. Atomic Physics 

We employ the so-called Mihalas-Hummer-Dappen 
(MHD) equation of state (EOS) (Hummer & Mihalas 
1988; Mihalas et al. 1988; Dappcn et al. 1988) which in- 
cludes explicit dissociation, ionization and excitation of 
all states of all species of all included elements. The only 
two molecular species included are H2 and H^, however. 
We have custom calculated tables of this EOS for the 
chemical mixture listed in Table 1, which includes 15 
elements as opposed to the previously published 6 ele- 
ment mixtures. As the dependent variables of the Navier 
Stoke's equations are g and e, we invert the EOS table to 
one in (lng,e). The entropy is evaluated by integration 
of the table 

S = S + |i(d £ -^dln^ , (2) 

where we first integrate along iso-chores of the table and 
then along the internal energy axis. The integration con- 
stant, So, is (arbitrarily) chosen so as to have zero-point 
at g = VTO x 10~ 7 gcm~ 3 and e = 2.5 x 10 12 ergcm~ 3 . 

The bound-free and free-free opacities, as well as scat- 
tering cross-sections are calculated using the MARCS 
package of codes (Gustafsson 1973), which is the engine 
behind the MARCS stellar atmosphere models (Gustafs- 
son et al. 1975, 2008). The data-tables, however, are 
updated as follows: 

We have included the absorption by bound-free 
(bf) (Broad & Reinhardt 1976; Wishart 1979), frcc-frcc 
(ff) (Bell & Bcrrington 1987), H+ bf+ff (Standi 1994), 
H^ ff (Bell 1980), and the photo-dissociation of OH and 
CH (Kurucz et al. 1987). The most important levels of 
He I, C I, N I, O I, Na I, Mg I, Mg II, Al I, Si I, Ca I, Ca II 
and Fe I were included as simple analytical fits as listed 
by Mathisen (1984). Bound-free absorption by O I and 
C II were adopted from the MARCS (Gustafsson et al. 
2008) implementation of OP opacities (Butler & Zeippen 
1990; Yan & Seaton 1987). Most of these changes were 
described in more detail by Trampedach (1997). 

Rayleigh scattering is included as fits to the results by 
Gavrila (1967) for H I, Langhoff et al. (1974) for He I and 
by Victor & Dalgarno (1969) for H2. The usual Thomson 
scattering by free electrons is also included. 

The line opacity is supplied by the opacity distribu- 
tion functions (ODFs) of Kurucz (1992a, b), which in- 
clude 58 million lines from the first nine ions of all ele- 



Table 1 

Logarithmic abundances normalized to A(H) = 12. 
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12.00 


14 


Si 
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He 
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Mg 


7.54 
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Ni 


6.21 
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Al 


6.43 









ments through Ni and the diatomic molecules H2, CO, 
CH, CN, C 2 , NH, OH, MgH, SiH, SiO and TiO. 

2.2. Chemical Abundances 

The Solar abundances we use, are listed in Table 1 as 
logarithmic abundances normalized to A(H) = 12.00. In 
our adaptation of atomic physics, the abundances enter 
in three independent places; In the EOS calculation, the 
bf- and ff-opacities and in the line opacities supplied by 
the ODFs. Since we have calculated our own tables of 
the MHD EOS and perform our own calculation of bf- 
and ff-opacities using the MARCS package, we are free 
to choose the abundances for those two parts - only re- 
stricted in the choice of elements by the availability of 
the necessary atomic physics. 

For the ODFs, however, we have less control over 
the assumed composition, and we have therefore largely 
adopted the abundances prescribed by the available 
ODFs at the time of starting the simulation grid. Kurucz 
(1992a, b) has made ODF tables available for the abun- 
dances by Anders & Grevesse (1989, AG89 hereafter) for 
a range of metallicities, but also for a He-free mixture 
and a few tables for lower Fe abundances. This enabled 
us to perform interpolations to the helioseismically de- 
termined Solar He abundance of A(He) = 10.92 (down 
from the classic value of A(Rc = 11.00) Basu & Antia 
(1995); Basu (1998) and the more modern Fe abundance 
of 7.50 (Grevesse & Sauval 1998; Asplund et al. 2009) 
(down from 7.67). These interpolations resulted in the 
abundances listed in Tabic 1. The differences between 
ours and the recent abundances based on solar convec- 
tion simulations, (Asplund et al. 2009, AGSS'09), are 
rather small and not systematic, except for the decreased 
C, N and O abundances. These three elements, however, 
have little effect on atmospheric structure, as they sup- 
ply only little line or continuum opacity, and affect the 
EOS minimally compared to the ionization of hydrogen 
(they have ionization energies similar to H). The only 
significant influence is through the molecules of H, C, 
N and O. Solar ATLAS9 models by Castelli & Kurucz 
(2003, only on-line) show very small differences in struc- 
ture, though, between GS'98 and AGS'05: The latter is 
20 K warmer at logTR, oss = —6.0 decreasing linearly to a 
few Kelvins around the photosphere and is about 30 K 
cooler below logr Ross ~ 1. AGSS'09 has C, N and O 
abundances that are closer to the "classic" GS'98 values, 
than does AGS'05, and we expect the above differences 
to be even smaller. 

For these reasons we did not feel compelled to recom- 
pute the, at that time, mostly finished grid of simulations 
to adopt the new abundances. 
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2.3. Radiative Transfer 

The radiative transfer in the atmosphere, has been 
described in detail by Nordlund (1982); Nordlund & 
Dravins (1990); Stein & Nordlund (2003). The for- 
ward solution is performed on long characteristics, and 
the original Feautrier (1964)-techniquc (See also Mihalas 
1978), has been modified to give more accurate solutions 
in the optical deep layers Nordlund (1982). The angle 
to the vertical of an outgoing characteristic (or ray) is 
/i = cos 9 and the azimuthal angle is ip. 

The angular integrations are performed with N v = 4 
(^-angles and = 2 /x-angles with /i-points and - 
weights determined by Radau quadrature (Radau 1880; 
Abramowitz & Stegun 1964). This method ensures the 
highest accuracy when one end-point (the vertical) is to 
be included. For = 2 this method gives 9 = 0° and 
70° with weights \ and |. The 0-angles are rotated 
for each time-step to avoid developing preferred direc- 
tions. To perform radiative transfer on slanted rays, the 
whole simulation cubes are rotated and tilted by rotating 
and shifting each horizontal plane sideways, exploiting 
the periodic boundaries. This is performed with bi-cubic 
spline interpolation. The number of grid-points in the in- 
clined cubes is the same as in the normal, vertical cube, 
and the computational cost is therefore also the same for 
the vertical as for each (/it, <j>) inclined calculation of the 
radiative transfer. 

The effects of spectral lines are included through the 
method of opacity binning. This consists of grouping 
wavelengths together based on their opacity strength, 
or more precisely, the Rosseland optical depth, tr oss , at 
which the monochromatic optical depth is t x = 1. Wc 
use four bins, with each bin spanning a decade in opac- 
ity strength. The source function used for each bin is 
simply the Planck function, B X (T), summed up over the 
wavelengths, Xj, in each bin, i, 

BiWi = ^BxjWXj , Wi = ^ w Aj . . (3) 

The weights, w Xj , for each wavelength is simply the wave- 
length width of each ODF, multiplied by the weight of 
each bin in the ODF. This method of opacity binning en- 
sures that the radiative cooling and heating in the pho- 
tosphere of strong lines is included as well as that of the 
continuum. The radiative heating per volume is 

Qrad=47T£ J K X (J X — S X )d\ (4) 

~A>KQK Koss '^ j Xi{Ji-Bi)Wi , (5) 

i 

where J x is the zero'th angular moment of the specific 
intensity, x . The general source function, S\, has been 
replaced by the strict LTE counterpart, B\. This means 
scattering has not been included (see Skartlien 2000; 
Hayek et al. 2010, for effects of scattering and a con- 
sistent formulation of it), except as absorption in the 
Rosseland mean (see below). The wavelength integration 
of Eq. (4) is approximated by the summing over bins, i, 
reducing the problem by about four orders of magnitude 
compared to a monochromatic calculation. The relative 
opacity, Xi = Ki/K^ oss , defines the bin- wise optical depth 
Ti = J QKY{ OSS Xidz. For each bin, Ji and Bi will diverge 



f° r T i ^5 1j gi vm g r i se to nrs t cooling in the photosphere, 
and then heating higher up in the atmosphere as Bi fol- 
lows the temperature and Jj converges to a constant. 
The fraction of heating/cooling contribution from each 
bin depends largely on the fraction of the source func- 
tion, Bi(j{ ~ l)wi, contained in that bin. 

The bin membership of each wavelength is determined 
from a ID, monochromatic radiative transfer calculation, 
performed on the average stratification of the simulation, 
averaged over a whole number of p modes (to yield a 
stable average). The bin membership is found as the 
Rosseland optical depth at which the given wavelength 
has unity optical depth 

i = int[log 10 t Ross (t x = 1)] , i G {0, 1, 2, 3} , (6) 

where 'hit' denotes the nearest integer. Wavelengths that 
go outside this range are assigned to the continuum bin, 
i = 0, or the strong-lines bin, i = 3, respectively For 
expediency we choose to simply scale the opacity of the 
continuum bin, such that Ki = 10* k . The proper opacity 
average within each bin, will need to converge to the dif- 
fusion approximation at depth, and the free-streaming 
approximation in the high atmosphere. We therefore 
choose a bridging 

k q = e- 2TR °^°Kj + (1 - e - 2T *»" ■°)kross,o , (7) 

between the Rosseland mean of the continuum bin and 
mean-intensity weighted opacity. These are defined as 

1 dB Xj I dB Xj 

«Ross,o- 2^ Kx . +ax . dT WX >/ ^ dT WX > ' 
j(i=o) Aj Aj / j(i=0) 

(8) 

where j(i = 0) is the set of wavelengths forming the 
continuum bin, according to Eq. (6). This Rosseland 
mean includes both the monochromatic absorption, k\, 
and monochromatic scattering, cr A , coefficients. 

The mean-intensity weighted opacity, kj, is summed 
over all wavelengths, but altered to suppress optically 
deep wavelengths with a factor e~ Tx ^ 2 




w Xj . (9) 



Note how scattering is not included in kj, except through 
the optical depth, dr x = q(k x + a x )dz, in the suppres- 
sion factor. The mean-intensity weighted opacity, kj, 
depends directly on the ID radiative transfer solution 
for the mean intensity, J x . To cover the whole opacity 
table, we extrapolate kj/kr oss o from the ID calibration 
stratification to the rest of the table. This is done along 
fits in log T and log g to iso-optical depth contours of the 
simulation. These tables are therefore specific to each 
individual simulation. 

In plane-parallel, gray radiative transfer calculations, 
the atmosphere converges to an iso-therm with height 
(King 1956, 0.81119 x T c s), as there is no heating or 
cooling above r ~ 0.1. In a line blanketed atmosphere, 
on the other hand, there are separate cooling peaks from 
lines of a range of strengths as outlined above, result- 
ing in atmospheres with decreasing temperatures with 
height and various features arising from the underlying 
opacities. 
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In deep layers, with continuum optical depth t± > 300 
for all points in the plane, we have added the radiative 
flux and its associated heating, as calculated in the dif- 
fusion approximation. 

2.3.1. Radiative contributions to the EOS 

The EOS tables already include the radiative contri- 
butions in the diffusion approximation, in particular for 
energy, e^ p = aT 4 /g, and pressure, p^ p = § T 4 . The 
radiation density constant is a — 8tt 5 k 4 / (15c 3 h 3 ) . We 
use the ID, monochromatic calibration to evaluate the 
proper expressions in the atmosphere 



_ 4tt J 
c g 



_deep 



J 



£ rad £j 



and 



Prad 



c 



K 

deep Ji 

^rad -q 



(10) 



(ii) 



where K is the second angular moment of the specific 

intensity. We therefore add £Vad P (ll — 1) to the pressure 
of the EOS table and equivalently to the internal energy. 
The J/B- and if/B-ratios are extrapolated from the ID 
average to the rest of the table, as described above. 

2.4. Relaxing the Simulations 

A simulation for a new choice of (T c g, logg)-parameters 
is started from a previous simulation with similar param- 
eters. The physical dimensions of the simulation box is 
scaled by the ratio of gravitational accelerations and the 
average entropy structure is changed to result in a new 
T e g based on all the previous simulations of the grid. The 
behavior of the entropy in the asymptotically deep inte- 
rior, with atmospheric parameters, is shown in Fig. 1. 
This asymptotic entropy is also what we feed into the 
simulations through the upflows at the bottom bound- 
ary, as confirmed from exponential fits to the horizon- 
tally averaged upflow entropy. The boundary affects the 
entropy by prematurely pulling it up to the asymptotic 
value, over the bottom 4-5 grid-points (0.3-0.5 pressure 
scale-heights). This boundary effect on entropy is small, 
though — only 0.4-1.5% of the atmospheric entropy jump. 

If we adjust only log g (with the associated scaling of 
the size of the box), but keep the entropy unchanged, the 
new simulation will end up along the adiabat of the orig- 
inal simulation and at the new logg. From Fig. 1 we see 
that those adiabats are diagonals in the plot. Many of 
the simulations lie along such adiabats, as this is the sim- 
plest and fastest way of starting a new simulation. The 
scaling of the box, should conceivably be accompanied 
by some scaling of the velocities. It turns out, however, 
that a factor of 10 2 change in g results in only a factor of 
1.5 change in vertical velocities (1.3 for horizontal veloc- 
ities). Keeping the fluxes consistent through the change, 
by not changing the velocities, seemed a better approach. 
These simulations will slump or expand, necessitating a 
new optimization of the vertical scale and extent. 

If T g needs to be adjusted away from the starting sim- 
ulations adiabat, more complicated adjustments must be 
invoked. First we shift the average entropy to the new 
■S'max and linearly stretch the average entropy stratifica- 
tion from the bottom to the atmospheric entropy mini- 
mum, to match the entropy jump. The expected jump 
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Figure 1. The asymptotic entropy (arbitrary zero-point, see be- 
low Eq. [2]), >Smax/[10 8 erg g _1 K _1 ] , of the deep convection zone 
as function of stellar atmospheric parameters. The T e ff-scale is 
logarithmic. The entropy is indicated with colors as shown on the 
color bar, and the location of the simulations arc shown with black 
asterisks, except for the solar simulation which is indicated with a 
0. For this figure only, we also added the simulation number from 
table 2. We have over-plotted tracks of stellar evolution computed 
with the MESA-code (Paxton et al. 2011), for masses as indicated 
along each track. The dashed part shows the pre-main-sequence 
contraction, and a and initial helium abundance, Yo, were deter- 
mined from a calibration to the present Sun. 

and S'max are found from inter- /extra-polations in Figs. 1 
and 4 between the previous simulations. We assume the 
simulations to be homologous on a gas pressure scale, 
Psc = Pgas /Pgas(peak in pturb), normalized at the loca- 
tion of the maximal pturb/ptot-ratio. The whole simula- 
tion cube is therefore adjusted adiabatically by the same 
pressure factor, and then adjusted iso-barically to the 
new entropy stratification. Our method does not rely 
on linearity of the EOS, but solves numerically for en- 
tropy along pressure contours. In both cases the changes, 
A In g and Ae, are found from the average stratification 
only, but applied to the whole cube. 

With these new pressures and densities, we scale the 
vertical velocities, u z ,to result in the projected peak 
Pturb/Ptot-ratio. We then adjust the amplitude of the 
internal energy fluctuations (keeping all the carefully ad- 
justed averages unchanged) in order to reproduce the 
target convective flux. We find a hydrostatic z-scale by 
inverting the equation of hydrostatic equilibrium 



dP 
~dz 



gg 



-ftot. bot 



<IL> 



(12) 



and integrating from the bottom and up. This z-scale 
will be rugged and not optimal for resolving the hydro- 
and thermo-dynamics. The last step is therefore to com- 
pute an optimized z-scale and interpolate the simulation 
cubes to this. This procedure results in simulations that 
are rather close to their (quasi-static) equilibrium state, 
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minimizing the relaxation time. 

No matter how carefully such a scaled simulation has 
been constructed, it will still not be in its new (quasi- 
static) equilibrium configuration, since we do not know 
a priori what that equilibrium is — that is the whole rea- 
son we need to perform these simulations, after all. The 
natural state is the lowest energy state, so departure from 
equilibrium means an excess of energy. This extra energy 
is quickly spent on exciting p modes (sound waves) and 
over a few sound crossing times, the energy gets chan- 
neled into 3-4 of the lowest order modes of the simulation 
box (radial, as well as non-radial) . We extract surplus en- 
ergy from the simulation by damping the radial p modes 
by adding a term 

"mode .,, (gu z ) 

— ,Wlth W mo dc = —pr- i ( 13 ) 

£damp \f?/ 

to the radial part of the momentum equation, Eq. (lb). 
For the damping time-scale, idamp, we use 1.3 times 
the period of the fundamental p mode, P\. When the 
center-of-mass velocity of the simulation no longer dis- 
plays simple sinusoidal modes, the damping time is in- 
creased smoothly (exponentially) until the damping is 
insignificant at which point it is turned off. The simula- 
tion is run without damping for another 2P\ for which we 
compute horizontal averages of g and T. This constitutes 
the first iteration of the relaxation process. 

If the stratification has changed from the previous iter- 
ation, the opacity binning is re-calculated with the new 
stratification and we start a new iteration of the relax- 
ation process (from Eq. (13). Otherwise we can start 
the production runs, from which the results in this pa- 
per are based. Each production run covers at least lOPi 
and all snapshots (about 100-300) of each simulation are 
included in the averaging producing our results. 

Some further requirements for a relaxed simulation are 
that the fluxes have reached a statistically steady state, 
that the total flux is statistically constant, that no quan- 
tities (especially at the bottom boundary) show any sig- 
nificant drift, that the radiative heating is statistically 
constant and that the time-steps are like-wise. All these 
precautions ensure that we have a homogeneous set of 
simulations and that we can trust in the differences be- 
tween the simulations and therefore can trust the depen- 
dencies of derived quantities with T c g and logg. 

2.5. The Grid of Simulations 

Our 37 simulations span most of the convective range 
of the zero age main sequence (ZAMS) from T c g — 
4 300 K to 6 900 K (K5-F3) and up to giants of log g = 2.2 
between T eS = 5 000K and 4 200 K (K2-K5), as listed in 
table 2. The location of the simulations in the atmo- 
spheric HR-diagram is also shown in Fig. 1. 

Our grid of simulations is not regular, for a couple of 
practical reasons. The most important one being that 
Toff is not an actual parameter of the simulations. We 
only have indirect control over T c s, by adjusting the en- 
tropy of the upflows entering through the bottom bound- 
ary. Requiring particular values of T c ff would add an- 
other outer loop of iterations of the relaxation process 
(Sect. 2.4), for monitoring and adjusting the resulting 
T c ff- Since this takes extra (human) time we have not re- 
quired a regular grid and pre-determined values of T e g. 



This is the reason that many of the simulations lie along 
adiabats of the deep entropy, as shown in Fig. 1. Those 
simulations result from changes to \ogg (with the asso- 
ciated scaling of the size of the box) only, without ad- 
justing the entropy structure. The entropy structure, of 
course, changes during the subsequent relaxation of the 
simulation, but the asymptotic deep value is fixed. 

Not insisting on a regular grid, also has the advantage 
of being able to have higher resolution around the main 
sequence, where it is most needed. This makes it possible 
to resolve features that would be missed in a coarser, 
regular grid. 

2.5.1. Interpolating in the Grid 

Interpolation in such an irregular (log T, log g)-grid is 
a bit more complicated compared to other grids of stellar 
atmospheres. Fortunately the triangulation of arbitrarily 
distributed points in a plane is not a new discipline. We 
use the suite of Fortran77 routines developed by Rcnka 
(1984), which performs both the triangulation and sub- 
sequent interpolations by piece-wise cubic 2D functions. 
The routines can also return the partial derivatives of 
the interpolation. The triangulation routines supplied in 
IDL 3 are also based on Rcnka (1984), facilitating visual- 
ization and interactive debugging. 

It makes no sense to interpolate individual snapshot 
of whole simulation cubes between different simulations, 
since they will not be homologous. We therefore only in- 
terpolate between temporally and horizontally averaged 
structures, or between derived quantities, e.g., limb dark- 
ening coefficients or the intensity contrast of granulation. 

2.5.2. Averaging Procedures 

The averaging of the simulations is carried out in two 
ways. First, a horizontal averaging mapped onto a (hori- 
zontally averaged) column density scale in order to filter 
out the main effect of the radial p-modes excited in the 
simulations. We call this pseudo Lagrangian averaging 
and denote it by (. . The reason for the "pseudo" is 
the use of a ID average column density scale, as opposed 
to the local column density. 

The second method we employ, is an optical depth av- 
eraging performed by averaging the quantities in ques- 
tion, over the undulating iso-r surfaces. This average we 
denote by (. . .) T and it has been calculated for both the 
Rosseland optical depth, tr OS s and the 5 000 A monochro- 
matic optical depth, T5ooo- 

What kind of average, to use, depends on the problem 
at hand, and should be guided by the structure of the 
equations involved (before approximations and simplifi- 
cations are made). 

3. GENERAL PROPERTIES OF SURFACE CONVECTION 

Deep convection is very nearly adiabatic (to parts per 
million), whether it occurs in the core of stars, or in the 
envelope. Here, no theory of convection is necessary for 
determining the structure of the star, as the stratifica- 
tion just follows an adiabat. This property makes deep 
convection very attractive for EOS research (Dappen & 
Gough 1984; Dappen & Nayfonov 2000), since the adia- 

3 IDL® is a registered trademark of ITT Visual Information 
Solutions, USA. 
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Table 2 

Fundamental parameters for the 37 simulations, and a few derived quantities, 
sim MK class log g p^j p^ p^ -p^ p^ 6 max 6 Jump 



1 


K3 


4681 ± 19 


2.200 


1248.30 


652.410 


3.478- 


16 


695 


283.010 


78.495 


7.9539 


6.7934 


0. 


26223 


0.18555 


2 


K2 


4962 ± 21 


2.200 


1800.00 


958.370 


4.598- 


27. 


616 


382.270 


170.020 


11.3540 


9.5480 





29639 


0.18758 


3 


K5 


4301 ± 17 


2.420 


667.81 


409.900 


1.724- 


12 


.100 


136.080 


32.745 


3.9781 


3.2772 


0. 


19130 


0.17758 


4 


K6 


4250 ± 11 


3.000 


165.00 


78.271 


0.457- 


- 1. 


.815 


31.953 


3.805 


1.6903 


1.6152 


0. 


14336 


0.16793 


5 


K3 


4665 ± 16 


3.000 


169.47 


115.880 


0.649 


- 2 


.564 


36.060 


6.236 


3.1968 


2.7581 


0. 


16881 


0.17371 


6 


Kl 


4994 ± 15 


2.930 


275.67 


132.840 


0.591 


- 3 


.752 


50.140 


9.476 


5.1703 


4.4027 


0. 


20124 


0.17806 


7 


G8 


5552 ± 17 


3.000 


328.29 


160.950 


0.787- 


- 4 


.101 


52.785 


17.077 


8.7944 


7.5617 


0. 


26233 


0.19568 


8 


K3 


4718 ± 15 


3.500 


52.18 


30.682 


0.176- 


- 


.751 


11.005 


1.203 


1.7158 


1.7188 


0. 


13479 


0.15879 


9 


K0 


5187 ± 17 


3.500 


53.59 


36.937 


0.196 


- 0. 


.891 


11.372 


1.826 


3.1948 


2.8700 


0. 


17238 


0.17183 


10 


K0 


5288 ± 20 


3.421 


89.00 


38.936 


0.205- 


- 0. 


.958 


15.609 


2.391 


3.9879 


3.5133 


0. 


17966 


0.17845 


11 


F9 


6105 ± 25 


3.500 


71.01 


52.645 


0.220- 


- 1. 


.416 


17.123 


5.964 


9.2519 


8.0430 


0. 


27125 


0.20493 


12 


K6 


4205 ± 8 


4.000 


14.83 


6.342 


0.039 


- 0. 


.138 


2.753 


0.116 


-0.3701 


0.5958 


0. 


07684 


0.11957 


13 


K4 


4494 ± 9 


4.000 


14.83 


6.420 


0.046 


- 


.113 


2.981 


0.165 


0.1071 


0.7660 


0. 


09391 


0.13550 


14 


K3 


4674 ± 8 


4.000 


14.83 


7.437 


0.034 


- 


.204 


2.945 


0.216 


0.4342 


0.9427 


0. 


10296 


0.13596 


15 


K2 


4986 ± 13 


4.000 


16.50 


9.026 


0.058- 


- 0. 


.196 


3.273 


0.299 


1.0546 


1.3418 


0. 


11742 


0.14480 


16 


G6 


5674 ± 16 


3.943 


19.33 


12.785 


0.069- 


- 0. 


.285 


4.069 


0.663 


3.1889 


2.9713 


0. 


17173 


0.17943 


17 


F9 


6137 ± 14 


4.040 


21.40 


9.326 


0.041- 


- 0. 


.300 


4.019 


0.747 


4.9301 


4.4363 


0. 


20478 


0.20117 


18 


F4 


6582 ± 26 


3.966 


24.28 


18.246 


0.092 


- 


.421 


6.034 


2.087 


9.2444 


8.1239 


0. 


26844 


0.21381 


19 


F4 


6617 ± 33 


4.000 


22.45 


16.877 


0.085 


- 0. 


.389 


5.368 


1.936 


9.2407 


8.1188 


0. 


25049 


0.21202 


20 


K4 


4604 ± 8 


4.300 


7.43 


3.625 


0.017- 


- 0. 


.101 


1.405 


0.076 


-0.2001 


0.6473 


0. 


08236 


0.12048 


21 


Kl 


4996 ± 17 


4.300 


7.43 


3.772 


0.019- 


- 0. 


.096 


1.515 


0.113 


0.4301 


0.9887 


0. 


10066 


0.13029 


22 


Kl 


5069 ± 11 


4.300 


8.27 


3.531 


0.023- 


- 


.074 


1.615 


0.115 


0.5441 


1.0562 


0. 


10300 


0.13118 


23 


K0 


5323 ± 16 


4.300 


8.27 


4.135 


0.023- 


- 0. 


.096 


1.614 


0.145 


1.0413 


1.3973 


0. 


11662 


0.14201 


24 


Gl 


5926 ± 18 


4.295 


9.32 


5.473 


0.024 


- 0. 


.139 


1.872 


0.259 


2.6796 


2.6358 


0. 


15869 


0.18234 


25 


F5 


6418 ± 26 


4.300 


11.76 


5.373 


0.030 


- 0. 


.128 


2.221 


0.406 


4.9238 


4.4645 


0. 


20351 


0.20745 


26 


F2 


6901 ± 29 


4.292 


11.45 


8.334 


0.037- 


- 0. 


.232 


2.760 


0.978 


9.2408 


8.1679 


0. 


27642 


0.20984 


27 


K4 


4500 ± 4 


4.500 


4.69 


2.108 


0.013- 


- 0. 


.051 


0.886 


0.031 


-0.6094 


0.4978 


0. 


06359 


0.11046 


28 


K3 


4813 ± 8 


4.500 


4.69 


2.026 


0.010- 


- 0. 


.054 


0.949 


0.049 


-0.2102 


0.6576 


0. 


07991 


0.11609 


29 


K0 


5232 ± 12 


4.500 


4.69 


2.344 


0.011- 


- 0. 


.062 


1.017 


0.071 


0.4277 


1.0217 


0. 


10218 


0.12981 


30 


G5 


5774 ± 17 


4.438 


6.02 


3.476 


0.020- 


- 


.082 


1.237 


0.140 


1.7112 


1.9149 


0. 


13747 


0.16345 


31 


F7 


6287 ± 15 


4.500 


5.36 


3.498 


0.020- 


- 0. 


.074 


1.178 


0.179 


3.1827 


3.0738 


0. 


16893 


0.19410 


32 


F4 


6569 ± 17 


4.450 


8.33 


3.704 


0.018- 


- 0. 


.097 


1.538 


0.288 


4.9211 


4.4866 


0. 


20431 


0.20658 


33 


Kl 


5021 ± 11 


4.550 


4.18 


2.055 


0.012- 


- 0. 


.050 


0.850 


0.052 


0.0039 


0.7731 


0. 


08862 


0.11827 


34 


G9 


5485 ± 14 


4.557 


4.17 


2.244 


0.013- 


- 


.050 


0.849 


0.073 


0.7452 


1.2347 


0. 


11311 


0.13862 


35 


Gl 


5905 ± 15 


4.550 


4.65 


2.547 


0.015- 


- 


.063 


0.933 


0.105 


1.7012 


1.9247 


0. 


13971 


0.16830 


36 


K6 


4185 ± 3 


4.740 


2.70 


1.016 


0.006- 


- 0. 


025 


0.488 


0.009 


-1.6403 


0.3562 


0. 


03898 


0.07347 


37 


K4 


4531 ± 10 


4.740 


2.70 


1.140 


0.007- 





.025 


0.539 


0.014 


-0.8481 


0.4077 


0. 


05513 


0.09950 



*This table is published in its entirety in the electronic edition of the Astrophysical Journal. 
a The simulation numbers are also shown in Fig. 1. 

b Thc spectral class is only approximate, as it depends on the luminosity. 

c Width of the simulation domain, w. 

d Dcpth of the simulation domain, d. 

e Min. and max. vertical grid-spacing. 

'The typical granulation size, Agran (See Sect. 4). 

g The convective expansion of the atmosphere, A (See Sect. 3.2). 

h The unit of entropy is 10 8 ergg -1 K _1 , and the zero-point is specified below Eq. (2). 



bat is completely determined by the EOS, and the matter 
is fully mixed, greatly simplifying the analysis. 

At the top of convective envelopes, the situation is 
much more complicated, not the least because it over- 
laps with the photospheric transition from optical deep 
(diffusion approximation) to optical thin (free-streaming 
approximation) and is accompanied by an abrupt change 
in density gradient. We outline the details below (also 
see Stein & Nordlund 1989; Nordlund & Stein 1991b, 
1997; Nordlund ct al. 2009). 

The upflows are close to isentropic, rising from the 
deep interior where convection is adiabatic. They also 
exhibit only little turbulence, since they are expanding 
along the exponential density gradient, smoothing out 
any fluctuations. The smoothing effect gets stronger as 
the density gradient and velocities increase towards the 
surface. This density gradient also forces the upflow to 



continuously overturn into the downdrafts, in order to 
conserve mass under the constraint of hydrostatic equi- 
librium. This overturning occurs with a scale-height of a 
mass mixing length, 

£ m = IdlnF^W/drf 1 , (14) 

where Ff m (r) is the mass flux, gu z , averaged over the 
upflows only (Trampedach & Stein 2011, TS11). This 
definition means a fraction of e _1 of the upflow will have 
overturned into the downdrafts over the height interval, 
£ m . In the deeper parts of our simulations (the same that 
were analyzed by TS11), below the supcradiabatic peak, 
the mass mixing length is proportional to the pressure 
scale height, l m = a m H p (not the density scale- height, 
H e , nor the distance to the top of the convection zone, 
A, as has been proposed in the past). We show this in 



Fig. 2 for three of our simulations that span our grid. 



a. 3 



J 



2: T, 
26: f, 
37: t, 



= 4962, logg = 2.20, a m =1.77±0.06 
:690B, logg=4.30, a m =1.72±0.05 
^4531, logg=4.74, a m =2.20±0.04 
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Figure 2. The mass mixing- length, £ m , in units of pressure scale- 
height, Hp, for three of our simulations (See table 2), as presented 
by Trampedach & Stein (2011). The abscissa is the logarithmic 
total pressure normalized by the pressure, p e fj, at the depth where 
(T) = T e g . The dashed, horizontal lines show the range of £ m /H p 
determined for each of the three simulations. The vertical dot- 
ted lines show the location of the top of the respective convection 
zones. The four asterisks on each curve, show the position of the 
temperature snapshots shown in Fig. 3. 

TS11 finds values of the deep a m between 1 .67 and 2.20. 
The sharp rise in a m at the bottom, is due to boundary 
effects. The mass flux is the quantity displaying the most 
extreme boundary effects, being further amplified by the 
differentiation, in this case. 

This simple, but powerful concept, of convective flows 
being dominated by continuous overturning of adiabatic 
upflows due to mass conservation of flows along a density 
gradient, is also predictive. At depths where radiative 
losses are negligible, pure advection of the entropy deficit 
from surface cooling into an isentropic interior, would 
result in a local entropy deficit, S(z) — 5 max , that scales 
with the density ratio, g(z = 0)/g(z). This is confirmed 
to high accuracy by the simulations (below z = 1 Mm for 
the solar case). 

There is a peak in £ m /H p about a pressure-decade be- 
low the top of the convection zone, before it drops off to 
about 1 above the photosphere. The top of the convec- 
tion zone is evaluated as the depth at which the convec- 
tive flux (Fconv = Fh + -Fkin, the sum of enthalpy and 
kinetic fluxes) is zero, and is indicated with the vertical 
dotted lines in Fig. 2. This point is well-defined since 
-Fconv is negative in the overshoot region above the con- 
vection zone. For comparison, the superadiabatic gradi- 
ent peaks around p/p c ff ~ 1, straddling the photosphere. 

The sharp drop in £ m from its peak value, reflects a 
change in morphology of the flows. As the mixing-length 
decreases faster than the pressure scale-height the same 
upflows have to overturn over a decreasing height-scale, 
forcing the whole circumference, and not just the ver- 
tices, of the upflows to be efficient downflows — this 
is what we recognize as granulation in the solar photo- 
sphere. This granulation pattern is formed in the front 
side of the mass mixing-length peak, as shown in Fig. 
3. Here we show the change of convective morphology 
with depth, for simulation 37 (sedate convection in a cool 
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Figure 3. Horizontal planes of temperature for arbitrary snap- 
shots of the three simulations shown in Fig. 2, and at the heights 
indicated with asterisks in that figure. The top row shows disk- 
center, white-light intensity. For each row the images are scaled 
so that contrasts can be directly compared, but the magnitude 
between simulations can not. 

MS dwarf), 26 (vigorous convection in a warm MS star) 
and 2 (vigorous convection in a red giant). The top row 
shows the emergent granulation pattern in disk-center, 
white-light intensity. The third and fourth rows show 
approximately the top and bottom of the granules and 
the second and bottom rows show convection six grid- 
points above and below that, respectively, to illustrate 
the transition away from granular convection. The dis- 
played quantity is temperature, but it is very similar in 
internal energy. Vertical velocity, or velocity divergence 
turned out to be less well suited for finding granules. 
We (loosely) define the bottom of the granules as where 
the cooler inter-granular lanes start to disconnect and 
and the downflows concentrate into knots connected by 
weaker lanes. 

In the layer at the top of the granulation (third row) 
we see all of the granules in the most sedately convec- 
tive simulation # 37. The most vigorous convection, as 
seen in # 26, is accompanied by a sharper photospheric 
temperature drop and we only see a few granules at the 
peak temperature (brightest in Fig. 3, third row) with 
the rest of the granules simply peaking in slightly deeper 
layers. The second row of images shows the morphology 
at the top of the convection zone, which is characterized 
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by hot lanes in the middle of inter-granular lanes, and 
some granules being cooler than these hot lanes at this 
height. These hot lanes are optically thin, however, and 
contribute little to the emergent picture of granulation, 
while the contribution by the granules is just fading, as 
only the hottest parts of the strongest granules have their 
photospheres at this height. 

The above described formation of granules is a very 
generic result. Whenever convection takes place in a 
stratified medium and the mass mixing length under- 
goes a sudden drop, faster than the pressure scale-height 
(typically by some form of cooling) a granulation pat- 
tern will emerge. It occurs in simulations with simple 
gray radiative transfer (Kim et al. 1995; Robinson et al. 
2003), in simulations with a completely artificial cool- 
ing layer and with impenetrable upper boundary (as in 
the interior convection simulations of, e.g., Miesch et al. 
2000; Brun et al. 2004), and it occurs in laboratory exper- 
iments with Rayleigh-Bernard convection (Koschmieder 
1967; Meyer et al. 1992; Busse & Clever 1998), where the 
mixing length drops to zero at the surface of the fluid. 

In real stellar convection, however, the contrast be- 
tween the upflows and the downdrafts is controlled by 
the details of the radiative cooling in the photosphere. 
This is dominated by the highly temperature sensitive 
H opacity (oc T 9 ). The isentropic upflows therefore 
only start cooling very close (spatially) to the unity op- 
tical depth. In the cooler downdrafts, on the other hand, 
unity optical depth occurs much deeper (on a spatial 
scale) (Georgobiani et al. 2003) and the cooling also oc- 
curs over a larger height range. The cooled downflows are 
also continuously mixed with the high-entropy plasma 
overturning from the upflows. As they move against the 
density gradient, the downflows are also turbulent and 
they contain the majority of vorticity — horizontal at the 
periphery of downdrafts and vertical in the centers. This 
inherent asymmetry between the upflows and the down- 
drafts makes it rather difficult to treat in analytical ID 
formulations. The asymmetry also results in density dif- 
ferences between upflows and downdrafts, which gives 
rise to a kinetic energy flux which is negative and in- 
creasing to about 20% of the total flux at the bottom 
of our simulation domains. This kinetic flux is not ac- 
counted for in MLT formulations, but together with the 
enthalpy flux makes up the convective flux in the simula- 
tions. There are also 5-13% contributions from acoustic 
fluxes inside the convection zone (largest for cool giants) 
and minor viscous fluxes that peak at less than 1%, just 
below the superadiabatic peak. The density difference 
in turn implies a difference in filling factor between the 
up- and downflows (assuming hydrostatic equilibrium as 
realized to a large extent by the simulations) . The filling 
factor is about 60% of the area being in the upflows in 
the deep parts, fairly constantly with height and atmo- 
spheric parameters. It changes abruptly at the top of 
the convection zone, to be a bit below 50% above the 
photosphere. 

3.1. Convective Efficiency 

The efficiency of convection is intimately connected to 
the convective velocities and the temperature and den- 
sity contrasts. A larger atmospheric entropy jump (a 
larger superadiabatic peak) means less efficient convec- 
tion, as the deviation from an adiabat implies energy 



losses on the way. The ideally efficient convection would 
not lose any energy before reaching the top of the con- 
vection zone. Comparing two cases with differing effi- 
ciency, but same total flux, the less efficient will necessar- 
ily have higher convective velocities in order to transport 
the same flux. These larger vertical velocities also means 
larger horizontal, overturning velocities between the up- 
flows and the downdrafts, sustaining larger pressure con- 
trast between the two. The entropy jump is set up by 
increased cooling around the photosphere, which drives a 
larger temperature contrast between the upflows and the 
downdrafts (and therefore also a larger emergent inten- 
sity contrast between granules and inter-granular lanes, 
as seen in Fig. 11). This entropy jump is shown in Fig. 4 
for our grid of simulations (also see Trampedach 2012). 




6000 5000 

T e „/m 

Figure 4. As Fig. 1, but showing the atmospheric entropy jump 
Sj U mp/[10 8 erg g - K _1 ], a measure of convective efficiency, as 
function of stellar atmospheric parameters. 

We clearly see that cool dwarfs have the most efficient 
convection and it gets progressively less efficient as the 
limit of convective envelopes is approached (a diagonal 
in the plot, above the high-T G ff/low-<7 limit of our grid). 
Beyond this limit no significant flux is transported by 
convection near the surface. Comparing with the asymp- 
totic adiabat shown in Fig. 1 (or the turbulent pressure 
ratio in Fig. 6), we see a very similar behavior, although 
the entropy jump falls off faster as the limit of convective 
envelopes is approached. In terms of MLT formulations, 
a large convective efficiency is accomplished with a large 
mixing length, which carries entropy further before mix- 
ing with the surroundings. 

Two stars that are otherwise similar (in effect, the same 
log<? and asymptotic entropy, S max ), but have different 
atmospheric opacities, will also have different convective 
efficiencies. The radiative loss of energy around the pho- 
tosphere, will happen in the same range of optical depth, 
for the two stars, but for the star with lower opacity 
this range will correspond to a larger range in (spatial) 
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depth, as dr = qk&z. With similar convective velocities, 
this would lead to more time spent loosing energy to the 
radiation field, and consequently increasing the atmo- 
spheric entropy jump and hence, decreasing the convec- 
tive efficiency. A lower efficiency leads to higher veloci- 
ties, which will limit the radiative losses above, providing 
a negative feed-back and enabling an equilibrium state. 
The net flux and therefore T e s will most likely be affected 
by such an opacity change. One reason for lower atmo- 
spheric opacities, is of course a lower metallicity, which 
would then be expected to result in lower convective ef- 
ficiency. We do not explore variation in metallicity with 
our grid of simulations, but a recent investigation of solar 
like Kepler stars by Bonaca et al. (2012), finds a decreas- 
ing a with decreasing metallicity, as argued above. This 
opacity effect on convection is, of course, also relevant 
for assessing effects of opacity updates on models or of 
how composition altering processes could affect stars and 
maybe explain observations. 

3.2. Seismology and Surface Convection 

The so-called surface effect (Christensen-Dalsgaard & 
Thompson 1997), is a systematic and persistent differ- 
ence between observed and model frequencies, which is 
independent of the degree, I, of the modes. This surface 
effect is caused by differences between the star and the 
model in the layers around the upper turning-point of 
the p modes. This is also the layers where the photo- 
spheric transition from optical deep (diffusion approxi- 
mation) to optical thin (free streaming), and the transi- 
tion at the top of the convection zone from convective to 
radiative transfer of the flux, occur. The latter results in 
the most inefficient convection in that star, with large ve- 
locities, large turbulent pressure, large temperature fluc- 
tuations and the most superadiabatic convection. The 
photospheric transition enhances the convective fluctua- 
tions by rapid cooling of the upflows. All these effects, 
and the strong coupling between radiation and convec- 
tion, means that a ID, plane-parallel, MLT atmosphere 
is an inadequate description. Most of the assumptions 
that MLT build on are violated in this region and we 
therefore expect, and have seen, large systematic differ- 
ences between observed and model p mode frequencies. 
Since the problems occur in a thin layer around the upper 
turning point, all modes feel it equally, and the frequency 
differences are therefore mainly a function of frequency. 
The usual procedure is therefore to adjust the model fre- 
quencies by a smoothly varying function of frequency. 
This procedure seems to work well for then Sun, and 
has also been applied successfully to other stars Kjeld- 
sen et al. (2008). 

The main reason for the surface effect is a convective 
expansion of the atmosphere, which enlarges the acoustic 
cavity of the modes. This expansion is shown in Fig. 5 
for the grid of simulations. It is evaluated as the aver- 
age outward shift of the (total) pressure stratification be- 
tween (lnptot)L of each simulation and their correspond- 
ing ID stellar envelope model, from their respective pho- 
tospheres and up. The ID models are based on the same 
EOS and opacities, and have been matched to the aver- 
aged simulations, thereby calibrating the mixing-length 
parameter, a (Trampedach et al. 2013). This ensures the 
3D and ID models are comparable, with the only major 
difference being the treatment of convection. 
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Figure 5. As Fig. 1, but showing the logarithmic total convective 
expansion of the atmosphere, A conv , as function of stellar atmo- 
spheric parameters. 

About half of this expansion is due to the turbulent 
pressure contributing to the hydrostatic support. The 
maximum of the turbulent- to total-pressure ratio is 
shown in Fig. 6, and is found to depend on atmospheric 
parameters in a similar way as the entropy of the asymp- 
totic adiabat in Fig. 1. Fig. 6 spans pressure ratios from 
3.90% for the coolest dwarf to 29.6% for the warmest gi- 
ant (closely followed by 27.6% for the hottest dwarf in 
our grid. From the different T c g and \ogg dependencies 
in Figs. 5 and 6 it is clear that it is not the pressure ratio 
itself that determines the atmospheric expansion. As we 
did for Eq. (12), we invert the equation of hydrostatic 
equilibrium and now integrate over just the turbulent 
component of the pressure, to get the associated expan- 
sion 

Aturb(z) = / rT — . (15) 

The total pressure is composed of, e.g., P = P r ad + -fgas + 
Ptuxbj and the gas pressure can be further decomposed 
into contributions by individual elements, ions, electrons 
or molecules. The amount of atmospheric expansion by 
any of those pressure components can be computed in the 
same way as that by the turbulent pressure, in Eq. (15) 
above. The integration is carried out from the interior 
to the top of the domain. About half of A tU rb (^quarter 
of the full A conv ) is actually supplied by the overshoot 
region above the convection zone, which in most ID MLT 
models has no velocity fields. 

This turbulent component provides approximately half 
of the full expansion shown in Fig. 5. The rest is due 
to the so-called convective back-warming (Trampedach 
2002, Paper I), which is caused by the large convective 
temperature fluctuations in the photosphere, combined 
with the high temperature sensitivity of the opacity there 
(k(H~) oc T 9 ). This effectively limits the radiative losses 
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Figure 6. As Fig. 1, but showing the maximum (with depth) of 
the Pturb/(^turb + F gas )-ratio, as function of stellar atmospheric 
parameters. 

from the hot granules, to a larger extent than would have 
been the case for a corresponding ID model with the av- 
erage temperature stratification. Because the tempera- 
ture fluctuations are well past the linear regime of the 
opacities, the effect is not symmetrically opposite in the 
cooler inter-granular lanes, and a net back-warming re- 
sults, akin to the effect of line-blanketing (Gustafsson 
et al. 2008, Sect. 6.1). The back-warming results in a 
larger (gas) pressure scale-height, i.e., an expansion of 
the atmosphere. 

Rosenthal et al. (1999) found the expansion of the so- 
lar atmosphere, relative to a ID MLT model, to be a 
major ingredient of the surface effect. The effective gra- 
dient of the turbulent pressure, 7i,turb (Stein & Nordlund 
1991), however, was found to be of similar importance. 
We show the inverse of this, dlng/dlnpturb for the solar 
simulation in Fig. 7. From the top panel of this plot, 
we see how (ptuib(t))h vs. (g(t))h displays a great deal 
of scatter, increasing with height above the photosphere. 
Despite this, a linear fit at each height is well-defined 
and gives a rather steep positive gradient above the pho- 
tosphere, going through vertical in the photosphere (the 
reason we show l/jx^uib in the bottom panel), and be- 
ing steep and negative in the interior, converging towards 
vertical (with minimal scatter). The shape of l/7i.tm-b 
is rather universal for our simulations, mainly varying in 
the amplitude of the various components of the curve. 

The appropriate form of the total 71 can be gleaned 
from the manner in which 71 enters the oscillation equa- 
tions. From Eq. (1) of Christensen-Dalsgaard (2008) we 
see that 71 only enters through terms similar to the first 
term on the right- hand-side of Eq. (16) 



1 dlnp _ 1 dlnpg 



1 din pturb 



7i jto tdlnr 71 dlnr ' 7i, tU rb dlnr ' 16 
and we propose the total 71 to be given by the left-hand- 
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Figure 7. Top panel: Black points show horizontally averaged 
(PturbM)L vs - {(?(*)} L f° r each snapshot of the solar simulation 
time-series (for every third depth-point). The green lines show 
linear fits to these points, at each depth. Bottom panel: Inverse 
gradients of the fits from the top panel, as function of depth (upper 
axis) and density (lower axis, in common with the top panel. 

side. Such a 71, tot requires a knowledge of the total pres- 
sure, p g -l-pturb- We can also define a total 71 that gives 
the response in the total pressure, based on just the gas 
pressure 



din; 



7l, tot 



dlnr 



/ 



1 din p, 



turb 



1 dlnpg 
71 dlnr 71, tm-b dlnr 



(17) 



We show this expression together with the pure gas 71 for 
our solar simulation, in Fig. 8. The reduced 71, adopted 
by Rosenthal et al. (1999) as a limiting case 



71, red = 71- 



if'g 



Pg + Pturb 



(18) 



is also shown in Fig. 8, and is rather different from our 
Eq. (17). Rosenthal et al. (1999) found that using the 
gas 71 for both pressures gave frequencies closer to the 
observed solar ones, than using Eq. (18). With the rather 
different behavior of Eq. (17), there is a potential for im- 
proved agreement with helioseismic observations, when 
using Eq. (17). 

This calculation is only a first step, since it assumed 
no phase-lag between density and turbulent pressure, and 
since we assumed 7i,turb to be independent of frequency. 
Future work with our grid of simulations will address 
these very important issue, which has grown more ur- 
gent with the advances in asteroseismology and the re- 
cent results by Mathur et al. (2012). They studied 22 
solar-like stars, targeted by NASA's Kepler mission, and 
found magnitudes of the surface effect that are not imme- 
diately reconcilable with Fig. 5, suggesting that 7i,turb- 
effects are significant. 
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Figure 8. The total 71, including both thermo- and hydro- 
dynamic effects, for the solar simulation of our grid. The peak 
hi 71, tot reaches a value of 5.46. 

4. VARIATION OF GRANULATION WITH STELLAR 
PARAMETERS 

In this section we present the variation of granulation 
size and intensity contrast with respect to the depen- 
dent variables of the grid, T e g and logg. In the case of 
our own star, the sizes of granules can be measured di- 
rectly from both ground- and space-based observations, 
and agree well with simulations. The intensity contrast, 
on the other hand, is less well constrained as observations 
are affected by both seeing, limited aperture, diffraction 
and scattering from various parts of the telescope, fo- 
cus issues, etc., that can easily halve the observed con- 
trast. Attempting to model these effects in Hinode ob- 
servations, Danilovic et al. (2008) found general agree- 
ment with their MuRAM simulations, with which ours 
agree fairly well. We find a 16.3% RMS fluctuation in 
white light granulation for our solar simulation, com- 
pared to their 14.4% RMS contrast in 6 300 A monochro- 
matic light. Beeck et al. (2012) shows that the 5 000 A 
contrast agrees between the simulations, although the 
detailed distribution and granulation patterns do differ 
some. 

It is worth noting that granulation in many ways are 
more constrained by unresolved spectral observations 
(Nordlund et al. 2009), as can also be carried out for 
stars. This is due to the fact that convection in late-type 
stars sets up correlations between temperature and ve- 
locity, which results in C-shaped bisectors of emergent 
spectral lines (Gray 1982; Basturk et al. 2011). The de- 
tailed line shapes of our solar and Procyon simulations 
agree well with observations (Asplund et al. 2000b; Al- 
lende Prieto et al. 2002). 

4.1. Intensity Distribution of Granulation 

We have calculated 2D spatial power spectra of the 
white light, disk center, surface intensity for each of the 
simulations, and averaged the spectra over time (at least 
10 periods of the fundamental p mode excited in each 
simulation). These spectra are very broad, as can be seen 
from Fig. 9, but the maxima, marked by filled circles and 
defining the typical size of granulation, is never-the-less 



well-defined. The range of each spectrum reflects the 
horizontal resolution at the lower end, and the horizon- 
tal extent at the upper end. The similar relative loca- 
tion of the peak in each spectrum means the granules are 
similarly well resolved. This granular size, A gran , is pre- 
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Figure 9. The time averaged spatial spectra of emergent disk- 
center intensity for each of the simulations of the grid. The thin 
lines are the raw spectra, and the thick lines are the spectra 
smoothed with a Gaussian kernel with FWHM of 11 wavcnum- 
ber points. We used short-dashed lines for giants (log g < 2.1), 
long-dashed lines for turn-off stars and solid lines for dwarfs 
(logg > 4.1). The solar simulation is shown in red. The abscissa, 
d = 2/(v / 7rA;), refers to diameters of circles with the same area as 
the 1/k 2 square, where k is the wavenumber corresponding to the 
Fourier transform. The maxima of the smoothed spectra are in- 
dicated with circles, and the corresponding granular Sizes, Agraji, 
are listed in table 2. The transformation of the maxima to the 
atmospheric HR-diagram, Fig. 10, is bijective, with giants to the 
right and cool dwarfs at the lower left. 

sented in Fig. 10. The granular size, A glan , grows almost 
inversely proportional with the surface gravity and in- 
creases only slightly with T e g. A least squares fit results 



log ? -f^ ~ (1.3210 ± 0.0038) x logT off 
|MmJ 

- (1.0970 ± 0.0003) x log g + 0.0306 ± 0.0359 

where g is in cgs-units and the la uncertainties of the 
linear regression is indicated. So the sensitivity to T e g 
is slightly larger than to logg, but since the stellar 
range of T e g is so much smaller than the range in grav- 
ity, the gravity dependence is effectively the most im- 
portant, as also seen in Fig. 10. It is often specu- 
lated that the size of granules scales with the pressure 
scale-height (Schwarzschild 1975) or the mixing length 
in the atmosphere, but we find such factors ranging from 
A glan /Hp ~ 9-13 and A graD /(aH p ) ~ 5-9, respectively, 
both generally increasing from dwarfs to giants. 
The time-averaged histograms of surface intensities has 
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Figure 10. The typical size of granulation (logarithmic scale) as 
function of atmospheric parameters. The granular size is found as 
the maxima of the smoothed spectra in Fig. 9. 

a bi-modal structure, with most of the surface being cov- 
ered in bright upflows with a narrower distribution of 
intensities than the less bright, inter-granular lanes. 

The intensity distributions of granules are narrower, 
due to the upflows being nearly isentropic (See Sect. 3). 
Their widths are determined by the radiative losses suf- 
fered by the over-turning upflows at the top of the con- 
vective envelope. If the top is located in optically deep 
layers, the radiative losses well be relatively minor and 
both the width of the intensity distribution of the gran- 
ules, and the intensity contrast to the inter-granular lanes 
are small — The granulation is largely hidden from our 
view. 

If the granulation reaches into optically thin layers, 
on the other hand, the reverse is the case and we see 
both a broader distribution and a larger contrast. This 
is often referred to as "naked granulation" , a term coined 
by Nordlund & Dravins (1990). 

The down-flows, on the other hand, are darker and 
have a much broader distribution. These intensity distri- 
butions are well-approximated by a simple double Gaus- 
sian 

n(I) = / ie -(C'-W3) 2 + he -((i-h)/hf (20) 

Fits to that expression are displayed in Fig. 12. Here 
we also notice that the simulations display less of a dip 
between the two components, compared to the double 
Gaussian fit. This third component, is due to the break- 
ing up of granules, which occurs in most of the snap- 
shots shown here. Some of the granules gets split apart 
by wedges of cool down-flows, others, mainly the largest 
granules, have cool spots developing in their centers, sud- 
denly succumbing to the negative buoyancy and turning 
into a down-draft. This occurs when the upflow in a gran- 
ule (proportional to the area of the granule) exceeds the 
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Figure 11. As Fig. 1, but showing the contrast of granulation, 
evaluated as the RMS scatter of the normalized intensity, 1/(1). 



capacity for horizontal outflow from the granule (propor- 
tional to the circumference), braking the upflow and thus 
increasing the cooling time. The developing cool spot, 
quickly connects to the nearest inter-granular lane. This 
phenomenon is also known as exploding granules (Namba 
& van Rijsbergen 1977; Spruit et al. 1990; Nordlund & 
Stein 1991a), although "exploding" seems a bit dramatic 
for all, but the largest of these events. (See, e.g., Rast 
1995, for comparisons of simulations and observations of 
exploding granules). 

The granulation in the snapshots of Fig. 12 looks rather 
different between the simulations, but this is a random 
occurrence, and the morphology looks different for an- 
other set of snapshots of the same simulations (cf. Fig. 
3 for different snapshots of simulations 2, 26 and 37). 
The over-all scale of the granulation is not apparent in 
Fig. 12 as the dimensions of each simulation is chosen 
so as to contain a similar number of granules (about 30 
granules), and the only property that varies markedly 
and systematically with atmospheric parameters, is the 
intensity contrast. We note that the distribution appears 
bimodal for only the higher gravity simulations. For the 
two giants of Fig. 12, however, the distributions are not 
symmetric and they are still best described by two sepa- 
rate Gaussians, as in Eq. (20). It is interesting that it is 
rather hard to pick out the intensity distribution from the 
picture. One might recognize sharper edges in the gran- 
ules of simulation #26, compared to that of #4, which 
gives rise to the separation of the two components in the 
distribution, but it is not immediately obvious that, e.g., 
the two giants are similar. We also note that the un- 
fitted component from splitting or "exploding" granules, 
increases with surface gravity, to the point of being ab- 
sent in the simulations of giants. This behavior is not 
yet understood. 

The contrast of the granulation, evaluated as the RMS 
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scatter around the mean is shown in Fig. 11. From this 
figure we see how the contrast behaves similarly to the 
convective efficiency, as measured by the atmospheric en- 
tropy jump, Fig. 4, with the highest contrast for the 
least efficient convection. The slope of iso-contrast lines 
change significantly from the cool to the warm side of the 
diagram, as opposed to the near-constant slopes of the 
entropy jump. The granulation contrast increases rather 
linearly across the contours, whereas the entropy jump 
increases close to exponentially, as function of logT g 
and \ogg. 

On the main sequence, for log g > 4, a linear fit to the 



intensity contrast gives 

/rms/[%] = (54.98 ± 1.86) x logT eff 

- (4.80 ± 0.53) x log g - 169.00 ± 0.73 , v 
and towards the giants, for log g < 3, 



(21) 



W/[%] = (20.81 ± 1.34) x logTeff 

- (1.11 ± 0.14) x logg - 55.46 ± 0.29 



(22) 



Notice the change of logT c ff- and log g-dependency be- 
tween dwarfs and giants, recognized as the change of 
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slope of the contours in Fig. 11. 

5. CONCLUSIONS 

We present an attempt at constructing a homogeneous, 
comprehensive, grid of 3D convective atmosphere simu- 
lations, with derived quantities for use by the general as- 
tronomical community. The simulations are constructed 
to be as realistic as possible, with realistic EOS, opac- 
ities and radiative transfer, in order to provide a solid 
foundation for interpreting observations. 

For all the subjects studied in the present paper, we 
found significant differences with respect to conventional 
ID stellar atmosphere models — differences that will have 
an impact on the interpretation of most stellar observa- 
tions. 

In Sect. 3 we gave an overview of the nature and mor- 
phology of convection, as observed in the simulations of 
our grid. We see fast, turbulent, entropy deficient down- 
drafts, driven by radiative cooling at the surface, plow- 
ing through an ambient, isentropic, laminar upflow. This 
is rather different from the picture normally associated 
with the MLT formulation, of distinct, warm, rising con- 
vective elements, existing for a mere mixing-length. We 
hope that the more realistic concept of convection, based 
on hydrodynamic simulations, will provide an improved 
framework for further theoretical developments in sim- 
plified, but realistic, descriptions of convection. 

Components of the seismic surface effect, a systematic 
frequency shift of observed p modes compared to ID stel- 
lar model predictions, were explored in Sect. 3.2. Part of 
the effect is due to enlarged acoustic cavities, caused by a 
convective expansion of the atmosphere (with respect to 
ID models). This expansion depends smoothly on atmo- 
spheric parameters, ranging from 8.7 km for the coolest 
dwarf in our grid, 140 km for the Sun, to 170 Mm for the 
warmest giant. We also attempt to derive an effective 
adiabatic gradient for the turbulent pressure, which is 
another important component of the surface effect. This 
7i,turb behaves rather differently from previous assump- 
tions about its properties. From the form of the adiabatic 
oscillation equations, we propose the effective 71 should 
be a harmonic mean, weighted by the pressure gradients, 
Eq. (16). 

Convection manifests itself as granulation at the sur- 
face, but has only been directly observed for the Sun. 
In Sect. 4, we find that the size depends smoothly on 
atmospheric parameters, but is not simply proportional 
to atmospheric pressure scale-height, as has often been 
assumed for a lack of observational constraints. 

We have shown that employing realistic 3D convec- 
tion simulations in the interpretation of stellar observa- 
tions will affect all stages of the analysis, and also allow 
new questions to be asked — questions that have been be- 
yond simplified, ID models of convection. In the near 
future our grid of simulations will be used for calibrating 
the MLT mixing-length, and evaluate the excitation and 
damping of p modes, among other applications. 
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